General instructions for R Markdown

Questions we are asking:

  1. Is increased parasitism associated with decreased body condition?
  • This is complicated because we are proposing an observational study and the causality that we are attempting to describe can go in both directions. Fish that have high parasitism could be skinner, but also skinner fish could be more susceptible to increased parasitism. We won’t be able to parse this out using our dataset, but it is still an interesting question worth answering. We just need to be careful to not make any absolute statements based on this (which is why we use the language ‘associated’ instead of ‘affected by’).
  • There is the element of time that we can incorporate into this as well, which is to say that we could look at body condition over time and see what trends exist. If there are any non-linear trends in body condition through time (meaning it doesn’t just go down or just go up or just stay the same), we can then incorporate the variable of parasite presence/burden to see if there is any correlation between the non-linear changes in body condition and parasite counts. This would basically look like on a graph, you have a big drop in body size followed a period of time later by an increase in parasite presence in these fish, you could say something like skinnier fish are associated with an increased parasite burden, since we would see that some event caused body size to decrease and as a result there was more parasites in the fish. But again we have to be really careful to not imply causation if all we have is correlation. Chelsea said she would be surprised if we noticed this effect, but it could be worth looking at anyway.
  1. Is the type of parasite important in determining whether parasitism is associated with decreased body condition?
  • Parasites exist in different parts of the body based on what they are: monogenes will always be on the gills, adult trematodes and adult cestodes will always be in the intestines, metacercariae can be most anywhere on the body, etc.
  • To answer this question, we would need to go through the data and identify what parasites are in their adult stage versus their larval stage. This would involve a bit of research into what the different (larval or adult) stages are present for each of the ones we have found. You would then be also able to look at the literature to see where that parasite would live within the fish at that stage in its life cycle, and then we could make commentaries about how different stages/types of parasites could be associated with changes to body condition. For example, if we had a ton of adult cestodes in the intestine and found that fish to have a very bad condition score, we could talk about how adult cestodes require more nutrients based on their big size, etc.

Methods

Rules for discriminating between adult and larval:

  • For any parasite we have identified to a specific taxonomic level (family, genus, etc.), we will look at that life cycle in the literature and ID what stage we have.
  • Anything labeled as META or CYST are considered larval, as metacercaria and cysts are by definition non-final adult forms of parasites.
  • Anything labeled as TREM, NEM, ACANTH, MONO, COPE, or CEST are considered adults of those groups, as we were able to ID them to larger taxonomic groups because they were not in cysts.
  • lists in the code below are the complete lists of which parasites are considered adult or larval

Fulton’s Condition Factor Equation:

Variable Definition
K Fulton’s Condition Factor (<1=good, >1=bad)
W weight of host fish (in grams)
SL standard length of fish (in millimeters)

\[ K = W/SL^3 *100 \]

Code for Analyses

Install packages

Load packages

library(ggeffects)
library(tidyverse)
library(ggpubr)
library(caret)
library(tibble)
library(readr)
library(dbplyr)
library(lubridate)
library(ggplot2)
library(glmmTMB)

Load datasets

Giving adult and larval IDs so we can look at life stage compared to parasite burden

Pimephales vigilax - add columns to identify species as an adult or larval stage and new column for total parasite count

Pvig.life.stage <- P.vigilax.data %>%
  mutate(adult = rowSums(select(., 
                      "acanth_bigb", 
                      "acanth_bkr",
                      "nem_cona",
                      "nem_cap",
                      "nem_dich",
                      "nem_larv",
                      "nem_trbk",
                      "nem_unk",
                      "cest_godz",
                      "cest_comp",
                      "cest_ntg",
                      "cope_cl",
                      "cope_grab",
                      "cest_pea",
                      "cest_l",
                      "cest_vit",
                      "cest_unk",
                      "mono_gyro",
                      "mono_dact",
                      "mono_lg",
                      "trem_rnda",
                      "trem_unk"))) %>%
  mutate(larval = rowSums(select(.,
                      "meta_hs",
                      "meta_unk", 
                      "myx_go",
                      "myx_geom",
                      "myxo_bc",
                      "myx_thel",
                      "myxo_sbad",
                      "myxo_unk",
                      "myxo_up",
                      "trem_met",
                      "trem_sss",
                      "trem_buc",
                      "trem_metaf",
                      "trem_metags",
                      "trem_ms",
                      "trem_nea"))) %>%
  mutate(adult = as.integer(adult)) %>% 
  mutate(larval = as.integer(larval)) %>%
  mutate(psite_count = rowSums(select(.,
                                      "adult", 
                                      "larval"))) %>%
  janitor::clean_names()

Ictalurus punctatus - add columns to identify species as an adult or larval stage and new column for total parasite count

Ipun.life.stage <- I.punctatus.data %>%
  mutate(larval = rowSums(select(., 
                      "nem_cyst",
                      "trem_meta_gg",
                      "trem_meta_unk",
                      "trem_lg",
                      "nem_cyst",
                      "myx_tail",
                      "myx_geom"))) %>%
  mutate(adult = rowSums(select(.,
                      "acan_ostr",
                      "cest_comp",
                      "cest_meg",
                      "leech_kut",
                      "nem_nmts",
                      "mono_ip",
                      "mono_unk",
                      "nem_phar",
                      "nem_pty",
                      "nem_sp",
                      "nem_taf",
                      "nem_unk",
                      "nmorph",
                      "trem_2",
                      "trem_allo",
                      "trem_ble",
                      "trem_crep",
                      "trem_f",
                      "trem_lip",
                      "trem_megict",
                      "trem_smile",
                      "trem_unk"))) %>%
  mutate(adult = as.integer(adult)) %>% 
  mutate(larval = as.integer(larval)) %>%
  mutate(psite_count = rowSums(select(.,
                                      "adult", 
                                      "larval"))) %>%
  janitor::clean_names()

Notropis athernoides - add columns to identify species as an adult or larval stage and new column for total parasite count

Nath.life.stage <- N.atherinoides.data %>%
  mutate(larval = rowSums(select(., 
                      "nem_cyst",
                      "trem_larv",
                      "trem_meta",
                      "myx_sp", 
                      "myx_round",
                      "trem_cyst"))) %>%#
  mutate(adult = rowSums(select(.,
                      "cest_comp",
                      "cest_hyd",
                      "cest_lvs",
                      "cest_tri",
                      "cope_a",
                      "mono_all",
                      "mono_unk",
                      "nem_bran",
                      "nem_cap",
                      "nem_spky",
                      "nem_tfs",
                      "nem_twas",
                      "nem_unkn",
                      "trem_gb",
                      "trem_l"))) %>%
  mutate(adult = as.integer(adult)) %>% 
  mutate(larval = as.integer(larval)) %>%
  mutate(psite_count = rowSums(select(.,
                                      "adult", 
                                      "larval"))) %>%
  janitor::clean_names()

Modify datasets to refect only parasites that are able to occur in larval and adult stages in the host.

  • we are doing this because in order to have a truly comparable dataset of parasites and their association with host body condition, we need to analyze parasites in a uniform way that does not skew the data one direction or another. For example, if you look at trematode parasites, both metacercaria and adult trematodes can live within a fish host. Now, compare that to monogenes who only occur as adults in the fish host. Because their larval stage will not occur in these fish hosts, by including them in our analyses we will be skewing the adult subset of parasites high. The monogenes have no larval counterpart in the fish, so they will always push the adult count higher. By removing monogenes from our analyses, we will have an even ground to compare larval and adult parasites. Therefore, parasites to NOT include in our analyses are:
  • monogenes
  • copepods
  • acanthocephalans
  • myxozoans

Mutate the data to remove the above-listed parasite groups from the analyses.

####
Pvig.life.stage <- P.vigilax.data %>%
  mutate(adult = rowSums(select(., 
                      "acanth_bigb", 
                      "acanth_bkr",
                      "nem_cona",
                      "nem_cap",
                      "nem_dich",
                      "nem_larv",
                      "nem_trbk",
                      "nem_unk",
                      "cest_godz",
                      "cest_comp",
                      "cest_ntg",
                      "cope_cl",
                      "cope_grab",
                      "cest_pea",
                      "cest_l",
                      "cest_vit",
                      "cest_unk",
                      "mono_gyro",
                      "mono_dact",
                      "mono_lg",
                      "trem_rnda",
                      "trem_unk"))) %>%
  mutate(larval = rowSums(select(.,
                      "meta_hs",
                      "meta_unk", 
                      "myx_go",
                      "myx_geom",
                      "myxo_bc",
                      "myx_thel",
                      "myxo_sbad",
                      "myxo_unk",
                      "myxo_up",
                      "trem_met",
                      "trem_sss",
                      "trem_buc",
                      "trem_metaf",
                      "trem_metags",
                      "trem_ms",
                      "trem_nea"))) %>%
  mutate(adult = as.integer(adult)) %>% 
  mutate(larval = as.integer(larval)) %>%
  mutate(psite_count = rowSums(select(.,
                                      "adult", 
                                      "larval"))) %>% 
  select(-matches("^(myx_|myxo_|acanth_|mono_|cope_)")) %>% 
  janitor::clean_names()

####
Ipun.life.stage <- I.punctatus.data %>%
  mutate(larval = rowSums(select(., 
                      "nem_cyst",
                      "trem_meta_gg",
                      "trem_meta_unk",
                      "trem_lg",
                      "nem_cyst",
                      "myx_tail",
                      "myx_geom"))) %>%
  mutate(adult = rowSums(select(.,
                      "acan_ostr",
                      "cest_comp",
                      "cest_meg",
                      "leech_kut",
                      "nem_nmts",
                      "mono_ip",
                      "mono_unk",
                      "nem_phar",
                      "nem_pty",
                      "nem_sp",
                      "nem_taf",
                      "nem_unk",
                      "nmorph",
                      "trem_2",
                      "trem_allo",
                      "trem_ble",
                      "trem_crep",
                      "trem_f",
                      "trem_lip",
                      "trem_megict",
                      "trem_smile",
                      "trem_unk"))) %>%
  mutate(adult = as.integer(adult)) %>% 
  mutate(larval = as.integer(larval)) %>%
  mutate(psite_count = rowSums(select(.,
                                      "adult", 
                                      "larval"))) %>% 
  select(-matches("^(myx_|acan_|mono_|leech_|nmorph)")) %>%
  janitor::clean_names()

#### 
Nath.life.stage <- N.atherinoides.data %>%
  mutate(larval = rowSums(select(., 
                      "nem_cyst",
                      "trem_larv",
                      "trem_meta",
                      "myx_sp", 
                      "myx_round",
                      "trem_cyst"))) %>%#
  mutate(adult = rowSums(select(.,
                      "cest_comp",
                      "cest_hyd",
                      "cest_lvs",
                      "cest_tri",
                      "cope_a",
                      "mono_all",
                      "mono_unk",
                      "nem_bran",
                      "nem_cap",
                      "nem_spky",
                      "nem_tfs",
                      "nem_twas",
                      "nem_unkn",
                      "trem_gb",
                      "trem_l"))) %>%
  mutate(adult = as.integer(adult)) %>% 
  mutate(larval = as.integer(larval)) %>%
  mutate(psite_count = rowSums(select(.,
                                      "adult", 
                                      "larval")))  %>% 
  select(-matches("^(myx_|mono_|cope_)")) %>%
  janitor::clean_names()

Calculating Fulton’s Condition Factor

Nath.life.stage <- Nath.life.stage %>%
  mutate(body.index.ath = ((Nath.life.stage$weight_mg/((Nath.life.stage$standard_length_mm)^3)))*(100)) %>%
  janitor::clean_names()

Pvig.life.stage <- Pvig.life.stage %>%
  mutate(body.index.vig = ((Pvig.life.stage$weight_mg/((Pvig.life.stage$standard_length_mm)^3)))*(100)) %>%
  janitor::clean_names()

Ipun.life.stage <- Ipun.life.stage %>%
  mutate(body.index.pun = ((Ipun.life.stage$weight_mg/((Ipun.life.stage$standard_length_mm)^3)))*(100)) %>%
  janitor::clean_names()

Combine datasets and scale data. Add condition factor back into the dataframes.

na_ip_pv_combined <- bind_rows(Nath.life.stage, Pvig.life.stage, Ipun.life.stage, .id = "dataset")
na_ip_pv_pp <- preProcess(na_ip_pv_combined, method = c("center", "scale"))
na_ip_pv_scaled <- predict(na_ip_pv_pp, na_ip_pv_combined)
notath_scaled <- na_ip_pv_scaled %>% filter(dataset == "1") %>%
  mutate(body.index.notath = ((Nath.life.stage$weight_mg/((Nath.life.stage$standard_length_mm)^3)))*(100)) %>%
  janitor::clean_names()
pimvig_scaled <- na_ip_pv_scaled %>% filter(dataset == "2") %>%
  mutate(body.index.pimvig = ((Pvig.life.stage$weight_mg/((Pvig.life.stage$standard_length_mm)^3)))*(100)) %>%
  janitor::clean_names()
ictpun_scaled <- na_ip_pv_scaled %>% filter(dataset == "3") %>%
  mutate(body.index.ictpun = ((Ipun.life.stage$weight_mg/((Ipun.life.stage$standard_length_mm)^3)))*(100)) %>%
  janitor::clean_names()

preliminary body condition index plots using original dataframe (with above-mentioned exclusions)

Plot life stage data using original dataframe (with above-mentioned exclusions)

# plot condition factor vs binary code for life stage (aka "stage-specific parasite stuff")
NA.bc.a <- ggplot(Nath.life.stage,
       aes(x = adult, y = body_index_ath, color = sex)) + xlab("Adult")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

NA.bc.l <- ggplot(Nath.life.stage,
       aes(x = larval, y = body_index_ath, color = sex)) + xlab("Larval")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

# plot condition factor vs binary code for life stage (aka "stage-specific parasite stuff")
IP.bc.a <- ggplot(Ipun.life.stage,
       aes(x = adult, y = body_index_pun, color = sex)) + xlab("Adult")  + geom_point() + 
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

IP.bc.l <- ggplot(Ipun.life.stage, 
       aes(x = larval, y = body_index_pun, color = sex)) + xlab("Larval")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

# plot condition factor vs binary code for life stage (aka "stage-specific parasite stuff")
PV.bc.a <- ggplot(Pvig.life.stage, 
       aes(x = adult, y = body_index_vig, color = sex)) + xlab("Adult")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

PV.bc.l <- ggplot(Pvig.life.stage,
       aes(x = larval, y = body_index_vig, color = sex)) + xlab("Larval")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

ggarrange(IP.bc.a, IP.bc.l, NA.bc.a, NA.bc.l, PV.bc.a, PV.bc.a, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C", "D", "E", "F"))

# more variability than the others

## double check sex, weights
## should we exclude these outliers

preliminary body condition index plots using scaled dataframe (with above-mentioned exclusions)

Plot life stage datausing scaled dataframe (with above-mentioned exclusions)

# plot condition factor vs binary code for life stage (aka "stage-specific parasite stuff")
NA.bc.a <- ggplot(notath_scaled,
       aes(x = adult, y = body_index_ath, color = sex)) + xlab("Adult")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

NA.bc.l <- ggplot(notath_scaled,
       aes(x = larval, y = body_index_ath, color = sex)) + xlab("Larval")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

# plot condition factor vs binary code for life stage (aka "stage-specific parasite stuff")
IP.bc.a <- ggplot(ictpun_scaled,
       aes(x = adult, y = body_index_pun, color = sex)) + xlab("Adult")  + geom_point() + 
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

IP.bc.l <- ggplot(ictpun_scaled, 
       aes(x = larval, y = body_index_pun, color = sex)) + xlab("Larval")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

# plot condition factor vs binary code for life stage (aka "stage-specific parasite stuff")
PV.bc.a <- ggplot(pimvig_scaled, 
       aes(x = adult, y = body_index_vig, color = sex)) + xlab("Adult")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

PV.bc.l <- ggplot(pimvig_scaled,
       aes(x = larval, y = body_index_vig, color = sex)) + xlab("Larval")  + geom_point() +
  ylab("Fulton's Condition Factor")  + theme_bw() + stat_smooth(method=lm)

ggarrange(IP.bc.a, IP.bc.l, NA.bc.a, NA.bc.l, PV.bc.a, PV.bc.a, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C", "D", "E", "F"))

# more variability than the others

## double check sex, weights
## should we exclude these outliers

Models

We want to run a general linear model to test the research questions: 1. Is body condition associated with increased total parasite burden? - is there a linear association between increasing body condition and increasing parasite burden? - does FCF increase as parasite burden increases? 2. Is body condition associated with infection with adult or larval stages of parasites? - is there a linear association between increasing body condition and increasing abundance of larval or adult parasites in the host? - does FCF increase as adult or larval parasite presence increases? - variables: - predictor = parasite_count, life_stage –> discrete - response = FCF –> continuous - discrete predictor = ANOVA - multiple/ continuous and discrete predictors = ANCOVA?

Test Model Formula for question 1

\[ K_i\sim\alpha+\beta_1count+\epsilon_i \]

Test Model Formula for question 2

\[ K_i\sim\alpha+\beta_1stage+\beta_2count+\epsilon_i \] ### Run linear model to test above questions

# Question 1 (lm)
summary(FCF_ath_lm <- lm(body_index_ath ~ psite_count, data = Nath.life.stage))
## 
## Call:
## lm(formula = body_index_ath ~ psite_count, data = Nath.life.stage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8708 -0.2325 -0.1553 -0.0608  7.2784 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.088027   0.085561  12.716   <2e-16 ***
## psite_count 0.001006   0.003980   0.253    0.801    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9823 on 169 degrees of freedom
##   (37 observations deleted due to missingness)
## Multiple R-squared:  0.0003782,  Adjusted R-squared:  -0.005537 
## F-statistic: 0.06395 on 1 and 169 DF,  p-value: 0.8007
plot(FCF_ath_lm)

summary(FCF_ath_lm <- lm(body_index_pun ~ psite_count, data = Ipun.life.stage))
## 
## Call:
## lm(formula = body_index_pun ~ psite_count, data = Ipun.life.stage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8625 -0.5370 -0.3616 -0.1666 13.7487 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.860079   0.261188   7.122 6.59e-10 ***
## psite_count -0.002641   0.005931  -0.445    0.657    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.14 on 72 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.002746,   Adjusted R-squared:  -0.01111 
## F-statistic: 0.1982 on 1 and 72 DF,  p-value: 0.6575
plot(FCF_ath_lm)

summary(FCF_ath_lm <- lm(body_index_vig ~ psite_count, data = Pvig.life.stage))
## 
## Call:
## lm(formula = body_index_vig ~ psite_count, data = Pvig.life.stage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8966 -0.4336 -0.2843 -0.0619 12.7554 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.997946   0.305819   6.533 3.51e-08 ***
## psite_count -0.003891   0.008464  -0.460    0.648    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.874 on 49 degrees of freedom
##   (146 observations deleted due to missingness)
## Multiple R-squared:  0.004294,   Adjusted R-squared:  -0.01603 
## F-statistic: 0.2113 on 1 and 49 DF,  p-value: 0.6478
plot(FCF_ath_lm)

# Question 1 (ANOVA)
summary(FCF_ath_anova <- aov(body_index_ath ~ psite_count, data = Nath.life.stage))
##              Df Sum Sq Mean Sq F value Pr(>F)
## psite_count   1   0.06  0.0617   0.064  0.801
## Residuals   169 163.07  0.9649               
## 37 observations deleted due to missingness
plot(FCF_ath_anova)

summary(FCF_ath_anova <- aov(body_index_pun ~ psite_count, data = Ipun.life.stage))
##             Df Sum Sq Mean Sq F value Pr(>F)
## psite_count  1    0.9   0.907   0.198  0.657
## Residuals   72  329.6   4.578               
## 15 observations deleted due to missingness
plot(FCF_ath_anova)

summary(FCF_ath_anova <- aov(body_index_vig ~ psite_count, data = Pvig.life.stage))
##             Df Sum Sq Mean Sq F value Pr(>F)
## psite_count  1   0.74   0.742   0.211  0.648
## Residuals   49 172.13   3.513               
## 146 observations deleted due to missingness
plot(FCF_ath_anova)

# Question 2 (lm)
summary(FCF_ath_al_lm <- lm(body_index_ath ~ larval * adult * psite_count, data = Nath.life.stage))
## 
## Call:
## lm(formula = body_index_ath ~ larval * adult * psite_count, data = Nath.life.stage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9419 -0.2419 -0.1293 -0.0078  7.0639 
## 
## Coefficients: (1 not defined because of singularities)
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.017e+00  1.095e-01   9.290   <2e-16 ***
## larval                    4.827e-03  1.537e-02   0.314    0.754    
## adult                     1.220e-02  2.597e-02   0.470    0.639    
## psite_count                      NA         NA      NA       NA    
## larval:adult              5.821e-03  3.952e-03   1.473    0.143    
## larval:psite_count       -3.455e-05  9.992e-05  -0.346    0.730    
## adult:psite_count        -5.134e-04  6.986e-04  -0.735    0.463    
## larval:adult:psite_count -6.880e-05  5.109e-05  -1.347    0.180    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9831 on 164 degrees of freedom
##   (37 observations deleted due to missingness)
## Multiple R-squared:  0.02828,    Adjusted R-squared:  -0.007275 
## F-statistic: 0.7954 on 6 and 164 DF,  p-value: 0.5748
plot(FCF_ath_al_lm)

summary(FCF_ath_al_lm <- lm(body_index_pun ~ larval * adult * psite_count, data = Ipun.life.stage))
## 
## Call:
## lm(formula = body_index_pun ~ larval * adult * psite_count, data = Ipun.life.stage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9417 -0.5111 -0.3990 -0.1282 13.6316 
## 
## Coefficients: (1 not defined because of singularities)
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.9771812  0.3748861   5.274 1.54e-06 ***
## larval                   -0.2912635  0.5894503  -0.494    0.623    
## adult                    -0.0220916  0.0757006  -0.292    0.771    
## psite_count                      NA         NA      NA       NA    
## larval:adult              0.0190266  0.0422251   0.451    0.654    
## larval:psite_count        0.0040358  0.0079439   0.508    0.613    
## adult:psite_count         0.0002453  0.0016887   0.145    0.885    
## larval:adult:psite_count -0.0004097  0.0008118  -0.505    0.615    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.209 on 67 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.01098,    Adjusted R-squared:  -0.07758 
## F-statistic: 0.124 on 6 and 67 DF,  p-value: 0.9931
plot(FCF_ath_al_lm)

summary(FCF_ath_al_lm <- lm(body_index_vig ~ larval * adult * psite_count, data = Pvig.life.stage))
## 
## Call:
## lm(formula = body_index_vig ~ larval * adult * psite_count, data = Pvig.life.stage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9284 -0.4533 -0.2908 -0.0267 12.6762 
## 
## Coefficients: (1 not defined because of singularities)
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               2.104e+00  4.792e-01   4.390 7.01e-05 ***
## larval                   -2.445e-02  5.497e-02  -0.445    0.659    
## adult                     4.899e-02  3.044e-01   0.161    0.873    
## psite_count                      NA         NA      NA       NA    
## larval:adult              1.423e-02  4.344e-02   0.327    0.745    
## larval:psite_count        1.495e-04  3.699e-04   0.404    0.688    
## adult:psite_count        -1.132e-02  3.694e-02  -0.306    0.761    
## larval:adult:psite_count -2.054e-05  9.178e-05  -0.224    0.824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.972 on 44 degrees of freedom
##   (146 observations deleted due to missingness)
## Multiple R-squared:  0.01026,    Adjusted R-squared:  -0.1247 
## F-statistic: 0.07602 on 6 and 44 DF,  p-value: 0.9981
plot(FCF_ath_al_lm)

# Question 2 (ANOVA)
summary(FCF_ath_al_anova <- aov(body_index_ath ~ larval * adult * psite_count, data = Nath.life.stage))
##                           Df Sum Sq Mean Sq F value Pr(>F)
## larval                     1   0.01  0.0079   0.008  0.928
## adult                      1   0.15  0.1524   0.158  0.692
## larval:adult               1   0.01  0.0100   0.010  0.919
## larval:psite_count         1   0.80  0.8010   0.829  0.364
## adult:psite_count          1   1.89  1.8884   1.954  0.164
## larval:adult:psite_count   1   1.75  1.7530   1.814  0.180
## Residuals                164 158.52  0.9666               
## 37 observations deleted due to missingness
plot(FCF_ath_al_anova)

summary(FCF_ath_al_anova <- aov(body_index_pun ~ larval * adult * psite_count, data = Ipun.life.stage))
##                          Df Sum Sq Mean Sq F value Pr(>F)
## larval                    1    0.5   0.532   0.109  0.742
## adult                     1    1.1   1.054   0.216  0.644
## larval:adult              1    0.5   0.470   0.096  0.757
## larval:psite_count        1    0.2   0.245   0.050  0.823
## adult:psite_count         1    0.1   0.086   0.018  0.895
## larval:adult:psite_count  1    1.2   1.243   0.255  0.615
## Residuals                67  326.9   4.879               
## 15 observations deleted due to missingness
plot(FCF_ath_al_anova)

summary(FCF_ath_al_anova <- aov(body_index_vig ~ larval * adult * psite_count, data = Pvig.life.stage))
##                          Df Sum Sq Mean Sq F value Pr(>F)
## larval                    1   0.67   0.669   0.172  0.680
## adult                     1   0.31   0.305   0.079  0.781
## larval:adult              1   0.01   0.007   0.002  0.967
## larval:psite_count        1   0.37   0.375   0.096  0.758
## adult:psite_count         1   0.22   0.223   0.057  0.812
## larval:adult:psite_count  1   0.19   0.195   0.050  0.824
## Residuals                44 171.10   3.889               
## 146 observations deleted due to missingness
plot(FCF_ath_al_anova)